To quantitatively examine the efficacy of vegetation restoration in drylands globally.
Study-level viz to document patterns in exclusions primarily and the relatie frequenices, at the study level, of major categories of evidence.
#study data####
library(tidyverse)
studies <- read_csv("data/studies.csv")
studies
## # A tibble: 278 x 18
## ID title technique data region exclude rationale observations
## <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 152 Shor… seeding,… expe… Africa no <NA> <NA>
## 2 180 Rest… chemical… App.… Africa no <NA> <NA>
## 3 229 Infl… soil see… fiel… Africa no <NA> <NA>
## 4 230 Acti… planting fiel… Africa no <NA> <NA>
## 5 255 The … grazing … fiel… Africa no <NA> <NA>
## 6 262 Reve… seeding,… eper… Africa no <NA> <NA>
## 7 263 The … phytogen… fiel… Africa no <NA> <NA>
## 8 264 Eval… seeding,… fiel… Africa no <NA> <NA>
## 9 271 Patc… natural … fiel… Africa no <NA> <NA>
## 10 4 Fact… natural … App.… Africa no <NA> <NA>
## # ... with 268 more rows, and 10 more variables: disturbance <chr>,
## # system <chr>, goal <chr>, intervention <chr>, paradigm <chr>,
## # grazing <chr>, hypothesis <chr>, soil <chr>, benchmark <chr>,
## # notes <chr>
#quick look at rationale needed
exclusions <- studies %>%
filter(exclude == "yes")
#quick look at studies with paradigms
evidence <- studies %>%
filter(exclude == "no")
#library(skimr)
#skim(evidence)
#study-level viz#####
#exclusions
#ggplot(exclusions, aes(rationale, fill = region)) +
# geom_bar() +
# coord_flip() +
# labs(x = "rational for exclusion", y = "frequency") +
# scale_fill_brewer(palette = "Paired")
ggplot(evidence, aes(disturbance, fill = paradigm)) +
geom_bar(na.rm = TRUE) +
coord_flip() +
scale_fill_brewer(palette = "Paired") +
labs(y = "frequency")
#ggplot(evidence, aes(region, fill = paradigm)) +
# geom_bar(na.rm = TRUE) +
# coord_flip() +
# scale_fill_brewer(palette = "Paired") +
# labs(y = "frequency")
#ggplot(evidence, aes(data, fill = paradigm)) +
# geom_bar(na.rm = TRUE) +
#coord_flip() +
#scale_fill_brewer(palette = "Paired") +
#labs(y = "frequency")
#ggplot(evidence, aes(system, fill = paradigm)) +
# geom_bar(na.rm = TRUE) +
# coord_flip() +
# scale_fill_brewer(palette = "Paired") +
# labs(y = "frequency")
#ggplot(evidence, aes(goal, fill = paradigm)) +
#geom_bar(na.rm = TRUE) +
#coord_flip() +
#scale_fill_brewer(palette = "Paired") +
#labs(x = "outcome", y = "frequency")
#step 1 models####
#paradigm
derived.evidence <- evidence %>%
group_by(technique, data, region, disturbance, goal, paradigm) %>% summarise(n = n())
#active-passive split
m <- glm(n~paradigm, family = poisson, derived.evidence)
anova(m, test="Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 167 9.9147
## paradigm 1 0.045115 166 9.8696 0.8318
#region
m1 <- glm(n~paradigm*region, family = poisson, derived.evidence)
#m1
#summary(m1)
anova(m1, test="Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 167 9.9147
## paradigm 1 0.045115 166 9.8696 0.8318
## region 6 0.301367 160 9.5682 0.9995
## paradigm:region 6 0.213627 154 9.3546 0.9998
#outcome
m2 <- glm(n~paradigm*goal, family = poisson, derived.evidence)
#m1
#summary(m1)
anova(m2, test="Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 167 9.9147
## paradigm 1 0.045115 166 9.8696 0.8318
## goal 6 0.240941 160 9.6287 0.9997
## paradigm:goal 4 0.301480 156 9.3272 0.9897
#even split between active and passive evidence by all key categories
A summary of sort process using PRISMA.
library(PRISMAstatement)
prisma(found = 1504,
found_other = 5,
no_dupes = 1039,
screened = 1039,
screen_exclusions = 861,
full_text = 178,
full_text_exclusions = 101,
qualitative = 77,
quantitative = 66,
width = 800, height = 800)
Check data and calculate necessary measures.
#all data
#data <- read_csv("data/data.csv")
#data <- data %>%
#mutate(lrr = log(mean.t/mean.c), rii = ((mean.t-mean.c)/(mean.t + mean.c)), var.es = ((sd.t^2/n.t*mean.t^2) + (sd.c^2/n.c*mean.c^2)))
#other effect size estimates
#library(compute.es)
#data <- data %>%
#mutate(d=mes(mean.t, mean.c, sd.t, sd.c, n.t, n.c, level = 95, #cer = 0.2, dig = 2, , id = ID, data = data))
#check metafor
#data from ag and grazing studies that examined restoration in drylands
data <- read_csv("data/data.csv")
mydata<- data %>%
filter(disturbance %in% c("agriculture","grazing")) %>%
filter(!notes %in% "couldnt extract data")
#write.csv(mydata, file="mydata.csv")
mydata <- mydata %>%
mutate(lrr = log(mean.t/mean.c), rii = ((mean.t-mean.c)/(mean.t + mean.c)), var.es = ((sd.t^2/(n.t*mean.t^2)) + (sd.c^2/(n.c*mean.c^2)))) #use only lrr now
#mydata<- mydata %>%
# mutate(aridity.range = cut(aridity.index, breaks=c(0,5,10,20,25,69.92), labels=c("hyperarid","arid", #"semiarid", "moderately arid", "slightly humid"))) #categorization of aridity.index values, according to de Martone
Explore summary level data of all data. Explore aggregation levels that support the most reasonable data structure and minimize non-independence issues.
#evidence map####
require(maps)
world<-map_data("world")
map<-ggplot() + geom_polygon(data=world, fill="gray50", aes(x=long, y=lat, group=group))
#map + geom_point(data=paperdata, aes(x=long, y=lat)) +
#labs(x = "longitude", y = "latitude") #render a literal map, i.e. evidence map, of where we study the niche in deserts globally
#add in levels and color code points on map####
#map of 178 articles
map + geom_point(data=data, aes(x=long, y=lat, color = paradigm)) +
scale_color_brewer(palette = "Paired") +
labs(x = "longitude", y = "latitude", color = "")
#map of selected articles, agriculture and grazing disturbances
map + geom_point(data=mydata, aes(x=long, y=lat, color = paradigm)) +
scale_color_brewer(palette = "Paired") +
labs(x = "longitude", y = "latitude", color = "")
#aggregation####
se <- function(x){
sd(x)/sqrt(length(x))
}
data.simple <- mydata %>%
group_by(study.ID, paradigm, technique, measure.success) %>%
summarise(n = n(), mean.lrr = mean(lrr), mean.rii = mean(rii), mean.var = mean(var.es))
main.data <- mydata %>%
group_by(study.ID, paradigm, intervention, outcome) %>%
summarise(n = n(), mean.lrr = mean(lrr), mean.rii = mean(rii), mean.var = mean(var.es))
#EDA data####
simple.data <- mydata %>% group_by(study.ID, paradigm, technique, measure.success) %>% summarise(mean.rii = mean(rii), error = se(rii))
simple.data <- na.omit(simple.data)
parad.data <- mydata %>% group_by(study.ID, paradigm) %>% summarise(mean.rii = mean(rii), error = se(rii))
parad.data <- na.omit(parad.data)
tech.data <- mydata %>% group_by(study.ID, technique) %>% summarise(mean.rii = mean(rii), error = se(rii))
tech.data <- na.omit(tech.data)
success.data <- mydata %>% group_by(study.ID, measure.success) %>% summarise(mean.rii = mean(rii), error = se(rii))
success.data <- na.omit(success.data)
#active
active <- mydata %>%
filter(paradigm == "active")
#viz for aggregation####
disturbance.data <- mydata %>%
group_by(study.ID,disturbance, paradigm) %>%
count()
disturbance.data2 <- disturbance.data %>%
group_by(disturbance,paradigm) %>%
count()
ggplot(na.omit(disturbance.data2), aes(disturbance,nn, fill=paradigm))+
geom_bar(stat = "identity") +
coord_flip(ylim=0:44) +
scale_fill_brewer(palette = "Blues")
intervention.data <- active %>%
group_by(study.ID,intervention, paradigm) %>%
count()
intervention.data2 <- intervention.data %>%
group_by(intervention,paradigm) %>%
count()
#ggplot(na.omit(intervention.data2), aes(intervention,nn, fill=paradigm)) +
#geom_bar(stat = "identity") +
#coord_flip(ylim=0:44) +
#scale_fill_brewer(palette = "Blues")
technique.data <- mydata %>%
group_by(study.ID,technique, paradigm) %>%
count()
technique.data2 <- technique.data %>%
group_by(technique,paradigm) %>%
count()
technique.data3 <- technique.data2 %>%
group_by(technique) %>%
count()
aridity<- mydata %>%
group_by(continent) %>% summarize(mean.aridity=mean(aridity.index))
aridity
## # A tibble: 6 x 2
## continent mean.aridity
## <chr> <dbl>
## 1 Africa 9.04
## 2 Asia 20.7
## 3 Europe 26.8
## 4 North America 19.6
## 5 Oceania 16.4
## 6 South America 22.9
#table(mydata$aridity.index)
ggplot(na.omit(data.simple), aes(technique, n, fill = paradigm)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_brewer(palette = "Paired")
ggplot(na.omit(data.simple), aes(measure.success, n, fill = paradigm)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_brewer(palette = "Paired")
#better
ggplot(main.data, aes(intervention, n, fill = paradigm)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_brewer(palette = "Paired") +
labs(fill = "")
ggplot(main.data, aes(outcome, n, fill = paradigm)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_brewer(palette = "Paired") +
labs(fill = "")
Exploratory data analyses to understand data and QA/QC using Rii.
Meta and conventional statistical models to explore relative efficacy.
Step 5. Synthesis stats
#p-value meta
#library(metap)
#all data for metas but cleaned for na's
#all data for meta cleaned
mdata.all <- mydata %>%
filter(!is.na(lrr)) %>%
filter(!is.na(var.es)) %>%
filter(!is.na(n.t)) %>%
filter(!is.na(p)) %>%
filter(!is.na(intervention)) %>%
filter(is.finite(lrr)) %>%
filter(!is.na(exp.length)) %>%
filter(!is.na(MAP)) %>%
filter(!is.na(aridity.index))
#after filtering NAs we apply models to 40 articles
summary(mdata.all$var.es)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00076 0.15430 0.01441 39.06266
#plot(density(mdata.all$var.es))
top_n(mdata.all, -10, var.es)
## # A tibble: 10 x 50
## study.ID ID publication.year data.year exp.year exp.length
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 111 111. 2016 2003 2003 132
## 2 111 111. 2016 2003 2003 132
## 3 111 111. 2016 2003 2003 132
## 4 111 111. 2016 2003 2003 132
## 5 111 111. 2016 2003 2003 132
## 6 111 111. 2016 2003 2003 132
## 7 111 111. 2016 2003 2003 132
## 8 111 111. 2016 2003 2003 132
## 9 111 111. 2016 2003 2003 132
## 10 111 111. 2016 2003 2003 132
## # ... with 44 more variables: disturbance <chr>, focus <chr>,
## # technique <chr>, intervention <chr>, paradigm <chr>, hypothesis <chr>,
## # pathway <chr>, plant.species <chr>, target.plant <chr>,
## # measure.success <chr>, outcome <chr>, measured.factor <chr>,
## # treatment <chr>, control <chr>, unit <chr>, Nsites <dbl>, n.t <dbl>,
## # n.c <dbl>, ntotalsamples <dbl>, mean.t <dbl>, mean.c <dbl>,
## # sd.t <dbl>, sd.c <dbl>, se.t <dbl>, se.c <dbl>, p <dbl>, df <dbl>,
## # measure.dispersion <chr>, lat <dbl>, long <dbl>, continent <chr>,
## # country <chr>, ecosystem <chr>, elevation <dbl>, MAP <dbl>, MAT <dbl>,
## # aridity.index <dbl>, `potential.evaporation (mm)` <dbl>,
## # grazing <chr>, soil <chr>, notes <chr>, lrr <dbl>, rii <dbl>,
## # var.es <dbl>
top_n(mdata.all, 5, var.es)
## # A tibble: 5 x 50
## study.ID ID publication.year data.year exp.year exp.length disturbance
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 11 11.7 2015 2007 1956 600 agriculture
## 2 11 11.7 2015 2007 1956 600 agriculture
## 3 11 11.7 2015 2007 1956 600 agriculture
## 4 11 11.7 2015 2007 1956 600 agriculture
## 5 147 147. 2012 2005 2005 36 agriculture
## # ... with 43 more variables: focus <chr>, technique <chr>,
## # intervention <chr>, paradigm <chr>, hypothesis <chr>, pathway <chr>,
## # plant.species <chr>, target.plant <chr>, measure.success <chr>,
## # outcome <chr>, measured.factor <chr>, treatment <chr>, control <chr>,
## # unit <chr>, Nsites <dbl>, n.t <dbl>, n.c <dbl>, ntotalsamples <dbl>,
## # mean.t <dbl>, mean.c <dbl>, sd.t <dbl>, sd.c <dbl>, se.t <dbl>,
## # se.c <dbl>, p <dbl>, df <dbl>, measure.dispersion <chr>, lat <dbl>,
## # long <dbl>, continent <chr>, country <chr>, ecosystem <chr>,
## # elevation <dbl>, MAP <dbl>, MAT <dbl>, aridity.index <dbl>,
## # `potential.evaporation (mm)` <dbl>, grazing <chr>, soil <chr>,
## # notes <chr>, lrr <dbl>, rii <dbl>, var.es <dbl>
max(mdata.all$var.es)/min(mdata.all$var.es)
## [1] 4.800933e+14
table(mdata.all$study.ID)
##
## 11 13 51 64 66 68 69 77 87 88 95 98 104 109 111 112 115 119
## 84 26 24 45 24 12 10 3 48 117 15 2 128 6 249 7 21 24
## 121 135 142 147 152 154 161 167 170 179 206 207 210 223 231 239 247 251
## 108 30 12 32 7 12 4 18 3 12 63 15 27 15 7 4 198 8
## 255 256 263 264
## 5 6 20 9
#active only data for meta cleaned up
mdata <- mydata %>%
filter(paradigm == "active") %>%
filter(!is.na(lrr)) %>%
filter(!is.na(var.es)) %>%
filter(!is.na(n.t)) %>%
filter(!is.na(p)) %>%
filter(!is.na(intervention)) %>%
filter(is.finite(lrr)) %>%
filter(!is.na(exp.length))
#passive only data for meta cleaned up
mdatap <- mydata %>%
filter(paradigm == "passive") %>%
filter(!is.na(lrr)) %>%
filter(!is.na(var.es)) %>%
filter(!is.na(n.t)) %>%
filter(!is.na(p)) %>%
filter(!is.na(intervention)) %>%
filter(is.finite(lrr)) %>%
filter(!is.na(exp.length))
#aggregated active data for metas var estimated with central tendency
#note - could also bootstrap mean variance here instead of arithematic mean
simple.mdata <- mdata %>%
group_by(intervention) %>%
summarise(lrr = mean(lrr), var.es = mean(var.es), n = mean(n.t))
simple.mdata.2 <- mdata %>%
group_by(intervention, outcome) %>%
summarise(lrr = mean(lrr), var.es = mean(var.es), n = mean(n.t))
simple.mdata.var <- mdata %>%
group_by(intervention) %>%
summarise(lrr = mean(lrr), var.es = se(var.es), n = mean(n.t))
simple.mdata2.var <- mdata %>%
group_by(intervention, outcome) %>%
summarise(lrr = mean(lrr), var.es = se(var.es), n = mean(n.t))
#aggregated passive data for metas var estimated with central tendency
simple.mdatap <- mdatap %>%
group_by(intervention) %>%
summarise(lrr = mean(lrr), var.es = mean(var.es), n = mean(n.t))
simple.mdatap.2 <- mdatap %>%
group_by(intervention, outcome) %>%
summarise(lrr = mean(lrr), var.es = mean(var.es), n = mean(n.t))
simple.mdatap.var <- mdatap %>%
group_by(intervention) %>%
summarise(lrr = mean(lrr), var.es = se(var.es), n = mean(n.t))
simple.mdatap2.var <- mdatap %>%
group_by(intervention, outcome) %>%
summarise(lrr = mean(lrr), var.es = se(var.es), n = mean(n.t))
#metas with p-values####
#schweder(mdata$p)
#sumz(p, data = mdata)
#mdata %>%
#split(.$intervention) %>%
#purrr::map(~sumz(.$p))
#sumlog(mdata$p)
#metas with meta package on effect size measures####
library(meta)
#all data non-aggregated
#m <- metagen(lrr, var.es, studlab = ID, byvar = intervention, data = mdata)
#summary(m)
#funnel(m)
#metabias(m)
#forest(m, sortvar = intervention)
#t-tests if different from 0
tmu <- function(x){t.test(x, mu = 0, paired = FALSE, var.equal=FALSE, conf.level = 0.95)}
mdata.all %>%
split(.$paradigm) %>%
purrr::map(~tmu(.$lrr)) #note this uses arithmetic means not estimated means from random effect models
## $active
##
## One Sample t-test
##
## data: x
## t = 7.6083, df = 1101, p-value = 5.943e-14
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.2069479 0.3507825
## sample estimates:
## mean of x
## 0.2788652
##
##
## $passive
##
## One Sample t-test
##
## data: x
## t = -7.5438, df = 357, p-value = 3.824e-13
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.4408271 -0.2585133
## sample estimates:
## mean of x
## -0.3496702
#Propose we only use random effect models because of the diversity of studies, differences in the methods and samples that may introduce heterogeneity
m1 <- metagen(lrr, var.es, studlab = ID, byvar = paradigm, comb.fixed=FALSE, data = mdata.all)
summary(m1) #active is positive and passive is negative so should not mix
## Number of studies combined: k = 1460
##
## 95%-CI z p-value
## Random effects model 0.0766 [0.0654; 0.0879] 13.38 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 184527270028.74 [184527270027.80; 184527270029.68]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 49679407227636228414242816.00 1459 0
##
## Results for subgroups (random effects model):
## k 95%-CI
## paradigm = active 1102 0.2184 [ 0.2055; 0.2314]
## paradigm = passive 358 -0.3413 [-0.3753; -0.3073]
## Q tau^2 I^2
## paradigm = active 49671693556322550581035008.00 0.0450 100.0%
## paradigm = passive 4152232320954931871744.00 0.1047 100.0%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 911.23 1 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
mb1<- metabias(m1)
mb1
#should do a t-test right now of paradigm is diff from 0
metareg(m1,aridity.index+exp.length)
##
## Mixed-Effects Model (k = 1460; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value): 0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 33284207035718112903168.00
## R^2 (amount of heterogeneity accounted for): 2.22%
##
## Test for Residual Heterogeneity:
## QE(df = 1457) = 48495089651041292387352576.0000, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 3535.8262, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.8886 0.0150 59.3593 <.0001 0.8593 0.9179 ***
## aridity.index -0.0339 0.0007 -51.0096 <.0001 -0.0352 -0.0326 ***
## exp.length -0.0007 0.0000 -26.2679 <.0001 -0.0008 -0.0006 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#interventions
m2 <- metagen(lrr, var.es, studlab = ID, byvar = intervention, comb.fixed=FALSE, subset = paradigm == "active", data = mdata.all)
summary(m2)
## Number of studies combined: k = 1102
##
## 95%-CI z p-value
## Random effects model 0.2184 [0.2055; 0.2314] 33.02 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 212403086959.62 [212403086958.53; 212403086960.71]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 49671693556322550581035008.00 1101 0
##
## Results for subgroups (random effects model):
## k 95%-CI
## intervention = vegetation 779 0.1845 [0.1694; 0.1996]
## intervention = soil 248 0.3128 [0.2990; 0.3265]
## intervention = water addition 75 0.6409 [0.5539; 0.7279]
## Q tau^2 I^2
## intervention = vegetation 48393423300974332080029696.00 0.0443 100.0%
## intervention = soil 97513761686148203151360.00 0.0111 100.0%
## intervention = water addition 31132.18 0.1047 99.8%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 226.58 2 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#metabias(m)
metareg(m2, ~ aridity.index + exp.length)
##
## Mixed-Effects Model (k = 1102; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value): 0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 44122812844877792935936.00
## R^2 (amount of heterogeneity accounted for): 2.22%
##
## Test for Residual Heterogeneity:
## QE(df = 1099) = 48490971316520689160159232.0000, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 1836.2421, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.5607 0.0193 29.1251 <.0001 0.5230 0.5985 ***
## aridity.index -0.0237 0.0008 -29.5763 <.0001 -0.0253 -0.0221 ***
## exp.length 0.0009 0.0000 19.0317 <.0001 0.0008 0.0010 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3 <- metagen(lrr, var.es, studlab = ID, byvar = intervention, subset = paradigm == "passive", comb.fixed=FALSE, data = mdata.all)
summary(m3)
## Number of studies combined: k = 358
##
## 95%-CI z p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 4152232320954931871744.00 357 0
##
## Results for subgroups (random effects model):
## k 95%-CI
## intervention = vegetation 125 0.2654 [ 0.2067; 0.3241]
## intervention = grazing exclusion 29 0.1351 [ 0.0270; 0.2431]
## intervention = soil 204 -0.7583 [-0.8196; -0.6970]
## Q tau^2 I^2
## intervention = vegetation 4152209524073903423488.00 0.1047 100.0%
## intervention = grazing exclusion 238316232.18 0.0881 100.0%
## intervention = soil 14453104123321616.00 0.1990 100.0%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 595.91 2 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#metabias(m)
metareg(m3, ~ aridity.index + exp.length)
##
## Mixed-Effects Model (k = 358; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.1047 (SE = 0.0876)
## tau (square root of estimated tau^2 value): 0.3236
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 11696393769174456320.00
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 355) = 4152219788056932122624.0000, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 815.4196, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3781 0.0565 6.6935 <.0001 0.2674 0.4888 ***
## aridity.index 0.0015 0.0029 0.5063 0.6127 -0.0042 0.0072
## exp.length -0.0020 0.0001 -23.0381 <.0001 -0.0021 -0.0018 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#trying models by continents
#m4 <- metagen(lrr, var.es, studlab = ID, byvar = continent, subset = paradigm == "active", comb.fixed=FALSE, data = #mdata.all)
#summary(m4)
#m5 <- metagen(lrr, var.es, studlab = ID, byvar = continent, subset = paradigm == "passive", comb.fixed=FALSE, data = #mdata.all)
#summary(m5)
#use random effects mean and var estimate to plot
#do t-tests here too
#outcomes
m6 <- metagen(lrr, var.es, studlab = ID, byvar = outcome, subset = paradigm == "active", comb.fixed=FALSE, data = mdata.all)
summary(m6)
## Number of studies combined: k = 1102
##
## 95%-CI z p-value
## Random effects model 0.2184 [0.2055; 0.2314] 33.02 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 212403086959.62 [212403086958.53; 212403086960.71]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 49671693556322550581035008.00 1101 0
##
## Results for subgroups (random effects model):
## k 95%-CI
## outcome = soil 249 0.2204 [ 0.1558; 0.2849]
## outcome = plants 305 0.5071 [ 0.4936; 0.5206]
## outcome = animals 24 -0.1152 [-0.1155; -0.1148]
## outcome = habitat 524 0.0621 [ 0.0437; 0.0804]
## Q tau^2 I^2
## outcome = soil 35077220764051.67 0.2656 100.0%
## outcome = plants 97513760543782884868096.00 0.0111 100.0%
## outcome = animals 541696.64 <0.0001 100.0%
## outcome = habitat 48393423303339003634253824.00 0.0443 100.0%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 8647.81 3 0
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#metabias(m)
metareg(m6, ~ aridity.index + exp.length)
##
## Mixed-Effects Model (k = 1102; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value): 0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 44122812844877792935936.00
## R^2 (amount of heterogeneity accounted for): 2.22%
##
## Test for Residual Heterogeneity:
## QE(df = 1099) = 48490971316520689160159232.0000, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 1836.2421, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.5607 0.0193 29.1251 <.0001 0.5230 0.5985 ***
## aridity.index -0.0237 0.0008 -29.5763 <.0001 -0.0253 -0.0221 ***
## exp.length 0.0009 0.0000 19.0317 <.0001 0.0008 0.0010 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m7 <- metagen(lrr, var.es, studlab = ID, byvar = outcome, subset = paradigm == "passive", comb.fixed=FALSE, data = mdata.all)
summary(m7)
## Number of studies combined: k = 358
##
## 95%-CI z p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 4152232320954931871744.00 357 0
##
## Results for subgroups (random effects model):
## k 95%-CI Q
## outcome = habitat 104 0.1605 [ 0.0964; 0.2246] 4152172019980636258304.00
## outcome = plants 50 0.4438 [ 0.0345; 0.8532] 33314950066906688.00
## outcome = soil 204 -0.7583 [-0.8196; -0.6970] 14453104123321616.00
## tau^2 I^2
## outcome = habitat 0.1047 100.0%
## outcome = plants 2.1620 100.0%
## outcome = soil 0.1990 100.0%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 425.72 2 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#metabias(m)
metareg(m7, ~ aridity.index + exp.length)
##
## Mixed-Effects Model (k = 358; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.1047 (SE = 0.0876)
## tau (square root of estimated tau^2 value): 0.3236
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 11696393769174456320.00
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 355) = 4152219788056932122624.0000, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 815.4196, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3781 0.0565 6.6935 <.0001 0.2674 0.4888 ***
## aridity.index 0.0015 0.0029 0.5063 0.6127 -0.0042 0.0072
## exp.length -0.0020 0.0001 -23.0381 <.0001 -0.0021 -0.0018 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#t-tests for outcomes diff from 0 but you can see using the 95% CI what is different
#check metafor for interactions?? in these big models or are we ok??
#brainstorm on how to explore?? techniques
#save but cut all this.
#m.study <- metagen(lrr, var.es, studlab = study.ID, byvar = intervention, data = mdata)
#summary(m.study)
#funnel(m)
#metabias(m)
#forest(m, sortvar = intervention)
#aggregated data
#m1 <- metagen(lrr, var.es, studlab = intervention, data = simple.mdata)
#summary(m1)
#funnel(m1)
#metabias(m1)
#forest(m, sortvar = intervention)
#different var estimate
#m2 <- metagen(lrr, var.es, studlab = intervention, data = simple.mdata.var)
#summary(m2)
#funnel(m2)
#metabias(m2)
#forest(m, sortvar = intervention)
#m3 <- metagen(lrr, var.es, studlab = intervention, byvar = outcome, data = simple.mdata.2)
#summary(m3)
#funnel(m)
#radial(m3)
#metabias(m2)
#forest(m, sortvar = intervention)
#m4 <- metagen(lrr, var.es, studlab = intervention, byvar = outcome, data = simple.mdata2.var)
#summary(m4)
#funnel(m)
#radial(m4)
#metabias(m2)
#forest(m, sortvar = intervention)
library(metafor)
#With metafor it is possible to fit metaregression models (e.g., Berkey et al. 1995; van Houwelingen et al. 2002), that #is, linear models that examine the infuence of one or more moderator variables on the outcomes (Viechtbauer, 2010).
# I´ve installed the 'devel' version of metafor (https://wviechtb.github.io/metafor/#installation)
#checking and comparing Lrr and sampling variances values obtained with escalc() function and by functions manually defined (lines 153-155)
#metadata<-escalc(measure="ROM",m1i=mean.t,m2i=mean.c,sd1i=sd.t,sd2i=sd.c,n1i=n.t,n2i=n.c, #data=mydata,var.names=c("LRR","LRR_var"),digits=4)
#metadata <- metadata %>%
# filter(!is.na(LRR)) %>%
# filter(!is.na(LRR_var)) %>%
# filter(!is.na(n.t)) %>%
# filter(!is.na(p)) %>%
# filter(!is.na(intervention)) %>%
# filter(is.finite(lrr)) %>%
# filter(!is.na(exp.length)) %>%
# filter(!is.na(aridity.index))
#summary(metadata$LRR_var)
#plot(density(metadata$var.es))
#top_n(metadata, -10, LRR_var)
#top_n(metadata, 5, LRR_var)
#write.csv(metadata, "metadata.csv")
#mdata.veg <- mydata %>%
# filter(paradigm == "active") %>%
#filter(intervention=="vegetation") %>% #filtering by intervention
#filter(!is.na(lrr)) %>%
#filter(!is.na(var.es)) %>%
# filter(!is.na(n.t)) %>%
# filter(!is.na(p)) %>%
# filter(!is.na(intervention)) %>%
# filter(is.finite(lrr)) %>%
# filter(!is.na(exp.length))
mdata.mean<- mdata.all %>%
group_by(study.ID, aridity.index, exp.length, paradigm, intervention, outcome) %>%
summarize(mean_lrr= mean(lrr),
mean_var.es= mean(var.es)) #collapsing data to means
#mdata.mean<- mdata.mean %>%
# filter(!is.infinite(mean_var.es)) #filter out outliers that are infinite
mod.1 <- rma(yi=lrr, vi=var.es, data = mdata.all)
summary(mod.1)
##
## Random-Effects Model (k = 1460; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -2057.9849 4115.9697 4119.9697 4130.5407 4119.9780
##
## tau^2 (estimated amount of total heterogeneity): 0.8959 (SE = 0.0341)
## tau (square root of estimated tau^2 value): 0.9465
## I^2 (total heterogeneity / total variability): 100.00%
## H^2 (total variability / sampling variability): 96407346870.76
##
## Test for Heterogeneity:
## Q(df = 1459) = 7665324384897.3271, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0899 0.0252 3.5678 0.0004 0.0405 0.1393 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(0.0899) #estimate
## [1] 1.094065
mod.2 <- rma(lrr, var.es, mods= ~paradigm+aridity.index+exp.length -1, data = mdata.all)
summary(mod.2)
##
## Mixed-Effects Model (k = 1460; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -1928.4827 3856.9655 3866.9655 3893.3827 3867.0069
##
## tau^2 (estimated amount of residual heterogeneity): 0.7501 (SE = 0.0287)
## tau (square root of estimated tau^2 value): 0.8661
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 80126342989.09
##
## Test for Residual Heterogeneity:
## QE(df = 1456) = 6538151352462.4238, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## QM(df = 4) = 290.7440, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## paradigmactive 0.9059 0.0601 15.0662 <.0001 0.7880 1.0237 ***
## paradigmpassive 0.5340 0.0924 5.7816 <.0001 0.3529 0.7150 ***
## aridity.index -0.0337 0.0027 -12.4170 <.0001 -0.0391 -0.0284 ***
## exp.length -0.0003 0.0001 -2.2753 0.0229 -0.0005 -0.0000 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Evaluating the influence of moderators on the heterogeneity
#interventions
mod.3 <- rma(lrr, var.es, slab= ID, mods= ~intervention+aridity.index+exp.length -1, data = mdata.all, subset = paradigm == "active")
summary(mod.3)
##
## Mixed-Effects Model (k = 1102; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -1470.7522 2941.5045 2953.5045 2983.5065 2953.5816
##
## tau^2 (estimated amount of residual heterogeneity): 0.7689 (SE = 0.0339)
## tau (square root of estimated tau^2 value): 0.8769
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 108731387818.35
##
## Test for Residual Heterogeneity:
## QE(df = 1097) = 6497632603139.7705, p-val < .0001
##
## Test of Moderators (coefficients 1:5):
## QM(df = 5) = 237.3672, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb
## interventionsoil 0.7767 0.0959 8.0997 <.0001 0.5887
## interventionvegetation 0.4676 0.0869 5.3814 <.0001 0.2973
## interventionwater addition 0.9680 0.1197 8.0847 <.0001 0.7333
## aridity.index -0.0258 0.0034 -7.5588 <.0001 -0.0324
## exp.length 0.0013 0.0002 5.9590 <.0001 0.0009
## ci.ub
## interventionsoil 0.9646 ***
## interventionvegetation 0.6379 ***
## interventionwater addition 1.2027 ***
## aridity.index -0.0191 ***
## exp.length 0.0017 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#forest(mod.3, slab= "study.ID")
mod.4 <- rma(lrr, var.es, slab= ID, mods= ~intervention+aridity.index+exp.length -1, data = mdata.all, subset = paradigm == "passive")
summary(mod.4) #aridity.index and exp.length non significants
##
## Mixed-Effects Model (k = 358; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -375.9568 751.9137 763.9137 787.1125 764.1565
##
## tau^2 (estimated amount of residual heterogeneity): 0.4459 (SE = 0.0350)
## tau (square root of estimated tau^2 value): 0.6677
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 498110093.57
##
## Test for Residual Heterogeneity:
## QE(df = 353) = 37647543864.7879, p-val < .0001
##
## Test of Moderators (coefficients 1:5):
## QM(df = 5) = 275.3915, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb
## interventiongrazing exclusion 0.2630 0.1818 1.4465 0.1480 -0.0934
## interventionsoil -0.1622 0.3593 -0.4515 0.6516 -0.8664
## interventionvegetation 0.3119 0.1185 2.6329 0.0085 0.0797
## aridity.index 0.0010 0.0064 0.1495 0.8812 -0.0116
## exp.length -0.0010 0.0007 -1.5800 0.1141 -0.0023
## ci.ub
## interventiongrazing exclusion 0.6194
## interventionsoil 0.5420
## interventionvegetation 0.5441 **
## aridity.index 0.0135
## exp.length 0.0003
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.4b <- rma(lrr, var.es, slab= ID, mods= ~intervention -1, data = mdata.all, subset = paradigm == "passive")
summary(mod.4b)
##
## Mixed-Effects Model (k = 358; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -378.4691 756.9381 764.9381 780.4266 765.0524
##
## tau^2 (estimated amount of residual heterogeneity): 0.4478 (SE = 0.0350)
## tau (square root of estimated tau^2 value): 0.6692
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 497656476.14
##
## Test for Residual Heterogeneity:
## QE(df = 355) = 37664401504.9221, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## QM(df = 3) = 271.4606, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb
## interventiongrazing exclusion 0.1358 0.1250 1.0865 0.2773 -0.1092
## interventionsoil -0.7540 0.0472 -15.9653 <.0001 -0.8465
## interventionvegetation 0.2479 0.0632 3.9228 <.0001 0.1240
## ci.ub
## interventiongrazing exclusion 0.3807
## interventionsoil -0.6614 ***
## interventionvegetation 0.3717 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#outcomes
mod.5 <- rma(lrr, var.es, slab= ID, mods= ~outcome+aridity.index+exp.length -1, data = mdata.all, subset = paradigm == "active")
summary(mod.5)
##
## Mixed-Effects Model (k = 1102; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -1456.9639 2913.9278 2927.9278 2962.9237 2928.0307
##
## tau^2 (estimated amount of residual heterogeneity): 0.7478 (SE = 0.0331)
## tau (square root of estimated tau^2 value): 0.8648
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 105846008972.89
##
## Test for Residual Heterogeneity:
## QE(df = 1096) = 6497561383370.6289, p-val < .0001
##
## Test of Moderators (coefficients 1:6):
## QM(df = 6) = 269.7943, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## outcomeanimals 0.1963 0.2040 0.9626 0.3358 -0.2034 0.5961
## outcomehabitat 0.0193 0.1158 0.1663 0.8679 -0.2077 0.2462
## outcomeplants 0.6239 0.0789 7.9103 <.0001 0.4693 0.7785 ***
## outcomesoil -0.1153 0.1500 -0.7683 0.4423 -0.4093 0.1788
## aridity.index -0.0087 0.0042 -2.0854 0.0370 -0.0169 -0.0005 *
## exp.length 0.0020 0.0003 7.5039 <.0001 0.0015 0.0026 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.6 <- rma(lrr, var.es, slab= ID, mods= ~outcome+aridity.index+exp.length -1, data = mdata.all, subset = paradigm == "passive")
summary(mod.6) #aridity.index non significant
##
## Mixed-Effects Model (k = 358; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -373.2789 746.5579 758.5579 781.7567 758.8006
##
## tau^2 (estimated amount of residual heterogeneity): 0.4420 (SE = 0.0347)
## tau (square root of estimated tau^2 value): 0.6648
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 493779723.64
##
## Test for Residual Heterogeneity:
## QE(df = 353) = 37647543969.8621, p-val < .0001
##
## Test of Moderators (coefficients 1:5):
## QM(df = 5) = 283.0688, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## outcomehabitat 0.2340 0.1221 1.9161 0.0553 -0.0054 0.4734 .
## outcomeplants 0.5128 0.1465 3.5002 0.0005 0.2256 0.7999 ***
## outcomesoil -0.0521 0.3463 -0.1503 0.8805 -0.7308 0.6267
## aridity.index 0.0004 0.0064 0.0633 0.9495 -0.0121 0.0129
## exp.length -0.0012 0.0006 -1.8689 0.0616 -0.0025 0.0001 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.6b <- rma(lrr, var.es, slab= ID, mods= ~outcome+exp.length -1, data = mdata.all, subset = paradigm == "passive")
summary(mod.6b)
##
## Mixed-Effects Model (k = 358; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -373.8439 747.6878 757.6878 777.0343 757.8602
##
## tau^2 (estimated amount of residual heterogeneity): 0.4407 (SE = 0.0346)
## tau (square root of estimated tau^2 value): 0.6638
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 490927793.15
##
## Test for Residual Heterogeneity:
## QE(df = 354) = 37647545438.6352, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## QM(df = 4) = 283.8777, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## outcomehabitat 0.2397 0.0839 2.8570 0.0043 0.0753 0.4042 **
## outcomeplants 0.5185 0.1142 4.5422 <.0001 0.2948 0.7422 ***
## outcomesoil -0.0518 0.3458 -0.1498 0.8809 -0.7296 0.6260
## exp.length -0.0012 0.0006 -2.0491 0.0405 -0.0023 -0.0001 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#forest(m1, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))
#forest(m2, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))
#forest(m3, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))
#forest(m4, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))
#t-tests for lrr diff from 0
#mdata %>%
#split(.$intervention) %>%
#purrr::map(~sumz(.$lrr))
#effect sizes plots
#need better ci estimates
#active plots
ggplot(simple.mdata, aes(intervention, lrr)) +
ylim(c(-2,2)) +
geom_point(position = position_dodge(width = 0.5)) +
labs(x = "", y = "lrr") +
coord_flip() +
geom_errorbar(aes(ymin=lrr-var.es, ymax=lrr+var.es), width=.05, position = position_dodge(width = 0.5)) +
geom_hline(yintercept = 0, colour="grey", linetype = "longdash")+
theme(axis.text.x=element_text(face="bold"),
axis.text.y=element_text(face="bold"),
axis.title=element_text(size=12,face="bold"),
strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold"))
ggplot(simple.mdata.2, aes(intervention, lrr, color = outcome)) +
ylim(c(-2,2)) +
geom_point(position = position_dodge(width = 0.5)) +
labs(x = "", y = "lrr", color = "") +
coord_flip() +
geom_errorbar(aes(ymin=lrr-var.es, ymax=lrr+var.es), width=.07, position = position_dodge(width = 0.5)) +
geom_hline(yintercept = 0, colour="grey", linetype = "longdash")+
theme(axis.text.x=element_text(face="bold"),
axis.text.y=element_text(face="bold"),
axis.title=element_text(size=12,face="bold"),
strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold"))
#passive plots
ggplot(simple.mdatap, aes(intervention, lrr)) +
ylim(c(-2,2)) +
geom_point(position = position_dodge(width = 0.5)) +
labs(x = "", y = "lrr") +
coord_flip() +
geom_errorbar(aes(ymin=lrr-var.es, ymax=lrr+var.es), width=.08, position = position_dodge(width = 0.5)) +
geom_hline(yintercept = 0, colour="grey", linetype = "longdash")+
theme(axis.text.x=element_text(face="bold"),
axis.text.y=element_text(face="bold"),
axis.title=element_text(size=12,face="bold"),
strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold"))
ggplot(simple.mdatap.2, aes(intervention, lrr, color = outcome)) +
ylim(c(-2,2)) +
geom_point(position = position_dodge(width = 0.5)) +
labs(x = "", y = "lrr", color = "") +
coord_flip() +
geom_errorbar(aes(ymin=lrr-var.es, ymax=lrr+var.es), width=.08, position = position_dodge(width = 0.5)) +
geom_hline(yintercept = 0, colour="grey", linetype = "longdash")+
theme(axis.text.x=element_text(face="bold"),
axis.text.y=element_text(face="bold"),
axis.title=element_text(size=12,face="bold"),
strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold"))
#ggplot(mdata, aes(lrr, color = intervention)) +
# geom_freqpoly(binwidth = 0.5, size = 2) +
# xlim(c(-10, 10)) +
# labs(x = "lrr", y = "frequency", color = "") +
# geom_vline(xintercept = 0, colour="grey", linetype = "longdash") +
# scale_color_brewer()
#ggplot(mdata, aes(lrr, fill = intervention)) +
# geom_dotplot(binwidth = 1) +
#xlim(c(-10, 10)) +
# labs(x = "lrr", y = "frequency", fill = "") +
# geom_vline(xintercept = 0, colour="grey", linetype = "longdash") +
# scale_fill_brewer()